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Abstract. 

We present Monte Carlo simulations of the extra galactic population of 
inspiralling double neutron stars, and estimate its contribution to the astrophysical 
gravitational wave background, in the frequency range of ground based interferometers, 
corresponding to the last thousand seconds before the last stable orbit when more than 
96 percent of the signal is released. We show that sources at redshift z > 0.5 contribute 
to a truly continuous background which may be detected by correlating third generation 
interferometers. 



1. Introduction 

Double neutron stars (DNSs) are very promising sources of gravitational waves for both 
ground based interferometers, which are sensitive to the last phase of the coalescence and 
the collapse, and the space detector LISA, which is expected to detect the continuous 
low frequency inspiral phase. In a previous work, we presented new estimates of the 
merging rate in the local Universe (including the contribution from ellipticals and taking 
into account the star formation history derived directly from observation) and discussed 
its consequences for the first generations of ground based interferometers. We predicted 
a detection every 148 and 125 years with Virgo and LIGO in their initial configuration, 
and up to 6 detections per year in their advanced configuration [D [2]- In addition 
to the emission from the nearest DNSs, it is expected that the superposition of a 
large number of unresolved sources produces a stochastic background. The background 
from the low frequency inspiral phase of various populations of compact binaries, which 
represents the main source of confusion noise for LISA, was studied intensively in the 
past decades (see for instance [31 IH El El E] for the galactic foregrounds, and [SI El UHl E] 
for the extra-galactic contribution). These predictions usually rely on binary evolution 
codes to estimate initial parameters such as eccentricity, mass ratio, orbital separation, 
which introduce large uncertainties due to the difficulty modelling mass loss and mass 
exchanges. In this work, we investigate the stochastic background produced by the extra- 
galactic population of DNSs, carrying on with the the previous estimates of [121 [13] ■ 
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We are interested in the last phase of the coalescence, up to the last stable orbit (LSO), 
at kHz frequencies, when the system is in circular orbit and when most of the GW 
energy is released. The article will be organized as follow: In section 1, we present 
a direct calculation of the spectral properties of the stochastic background; in section 
2, we discuss the detection of the background with different generations of terrestrial 
interferometers; in section 3, we present Monte Carlo simulations of the extra-galactic 
population of DNSs and its contribution to the GW background; and finally in section 
5, we summarize our results and discuss the interest of the Monte Carlo simulations for 
source modelling and data analysis, as well as possible improvements. 



2. The GW background 

The spectrum of the gravitational stochastic background is characterized by the 
dimensionless density parameter (or closure density) [HI [15] : 

= -p^ (1) 

where pgw is the gravitational energy density, Uo the frequency in the observer 
frame and pc = the critical energy density needed to close the Universe today. 
Throughout this paper, we assume a flat Einstein de Sitter 737 cosmology, with energy 
density of matterf2m = 0.3, energy density of vacuum = 0.7 and Hubble parameter 
-f^o = 70 km s~^ Mpc~^ [16], corresponding to the so-called concordant model derived 
from observations of distant type la supernovae [T7] and the power spectra of the cosmic 
microwave background fluctuations [18] . 

For a stochastic background of astrophysical origin 



^qw — q ^oFy (2j 

where the integrated flux (in erg cm^ Hz~^ s~^) at the observed frequency Vq is 
defined as: 

fj-^dz (3) 

The spectral properties of a single source located at redshift z are given by the 
fluence (in erg cm^ Hz~^) [191 120] : 

1 dEgyj 1 dEn 
i\ dUo A'nd\ 

where di = r{l + z) is the luminosity distance, r the proper distance, which depends 
on the adopted cosmological model, and u = Uoil+z) the frequency in the source frame. 
In the quadrupolar approximation, the spectral GW energy emitted by a binary system, 
which inspirals in a circular orbit is given by [2T1 122] : 

dEg^/dv = Kv~^''^ (5) 



47rrf? dv^ 47rrf? dv > ^ ^ 
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where 



K 



(Gtt)^/^ mim2 
3 (mi + 7712)-'^/^ 



(6) 



For double neutron stars with masses nii = m2 = 1.4 M©, one obtains K = 5.2 x 10^° 
erg Hz~^/^ and the gravitational frequency at the last stable orbit is assumed to be 



The merging of two neutron stars (NSs) occurs long after the formation of the 
progenitors and this delay must be taken into account when calculating the event rate. 
The progenitor formation rate per comoving volume is given on the time scale of the 
observer by: 



where R^{zf) is the cosmic star formation rate per comoving volume (SFR) at the time 
of formation Zf and is expressed in Mq Mpc~^ yr^^. In our calculations, we consider 
the recent model of [21], constrained by the Super Kamiokande limit on the electron 
antineutrino flux from past core-collapse supernovas up to Zmax = 6. The (1 + z) term 
in the denominator corrects the cosmic star formation rate by the time dilatation due to 
the cosmic expansion. The factor A is the mass fraction converted into the progenitors, 
assumed to be the same at all redshifts, and given by [2] as the product A = ^Nsh^NS, 
where (3ns is the fraction of binaries which remains bounded after the second supernova 
event, the fraction of massive binaries formed among all stars and Xns the mass 
fraction of NS progenitors derived from a modified Salpeter A IMF with minimal and 
maximal initial masses of 8 Mq and 40 Mq, as suggested by [23]. Numerically one 
obtains A = 3 x 10"^ M'^ 

The coalescence rate per comoving volume on the time scale of the observer results 
from the convolution between and the probability distribution of the coalescence 
time Tf. (or the delay between the formation of the DNS after the second supernova 
explosion and the coalescence), which depends on the orbital parameters (separation 
and eccentricity) and neutron star masses, at the time the DNS is formed, and is given 
by [2] as: 



ULso = 1.5 kHz m- 



R%Zf) = A 



R*{zf) 



(7) 



0.087 



with Tc e [0.2Myr - 20Gyr] 



(8) 



c 



In terms of the cosmic (lookback) time: 




(9) 



where 




(10) 



and where Zc is the redshift of coalescence, it writes: 



R^iz,) = i?°(T,) = I R%T, + n + T,)PrXrc)dT, 



(11) 
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Figure 1. comparison between the formation rate of DNS progenitors and the 
coalescence rate. 



where r?, (~ 10* yr) is the mean hfetime of the progenitors (or the time for the two 
massive stars to evolve into two neutron stars) |2|. The lower limit of the integral 
corresponds to the minimal coalescence time Tq = 0.2 Myr [2], while the upper limit is 
fixed by the maximal redshift of the adopted star formation rate and is given by the 
maximal formation time delayed according to the time to reach the coalescence, namely 

T(2;max) -n- Tc. FoY Z^^.^ = 6, T{z^^^) = Tmax = 12.5 Gy. 

As one can see in Fig. [T| the coalescence rate is shifted toward lower redshifts, with 
respect to the formation rate, reflecting the time delay between the formation of the 
progenitors and the coalescence event. 

The event rate per redshift interval in eq. [sj is thus given by multiplying R°{z) by 
the element of comoving volume: 



where 



'-f^K'^ (12) 

dz dz 

Combining the expressions above and replacing the constants by their usual values, 
one obtains: 

8.6 X , , dz (14) 

The upper limit of the integral, which depends on both the maximal emission 
frequency in the source frame Vmax and the maximal redshift of the model of star 
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formation history (-Zmax ~ 6), is given by: 



_ f 2:max if Z/o < 

™P otherwise ^ ^ 



I'd 



Consequently, the shape of the spectrum is characterized by a cutoff at the maximal 
emission frequency and a maximum at a frequency which depends on the shape of both 
the SFR and the spectral energy density. Before the maximum, ilg^, increases as u"^^^ 
(eq.[l5. 

Besides the spectral properties, it is important to study the nature of the 
background [25] • In the case of burst sources the integrated signal received at 2 = 
from sources up to redshift z, would show very different statistical behaviour whether 
the duty cycle 



D{z) = / f{l + z') — {z')dz' (16) 
Jo dz 

defined as the ratio, in the observer frame, of the typical duration of a single event 
f , to the average time interval between successive events, is smaller or larger than unity. 
When the number of sources is large enough for the time interval between events to be 
small compared to the duration of a single event {D » 1), the waveforms overlap to 
produce a continuous background. Due to the central limit theorem, such backgrounds 
obey the Gaussian statistic and are completely determined by their spectral properties. 
They could be detected by data analysis methods in the frequency domain such as the 
cross correlation statistic presented in the next section [26]. On the other hand, when 
the number of sources is small enough for the time interval between events to be long 
compared to the duration of a single event {D « 1), the sources are resolved and 
may be detected by data analysis techniques in the time domain (or the time frequency 
domain) such as match filtering [271 128]. An interesting intermediate case arises when 
the time interval between events is of the same order of the duration of a single event. 
These signals, which sound like crackling popcorn, are known as "popcorn noise". The 
waveforms may overlap but the statistic is not Gaussian anymore so that the amplitude 
on the detector at a given time is unpredictable. For such signals, data analysis strategies 
remain to be investigated [29], since the time dependence is important and data analysis 
techniques in the frequency domain, such as the cross correlation statistic, are not 
adapted. The critical redshifts at which the background becomes continuous, popcorn 
or shot noise will be fixed by the conditions D{zc) =1, 0.1 or 0.01 [25] . 

In our calculations, we considered the last ~ 1000 s before the last stable orbit, 
when more than 96% of the gravitational energy is released and when the signal range 
between 10 — 1500 Hz, in the frequency domain of ground based detectors [12]. At 
that time, the system has been circularized through GW emission (eq. |5]) and all the 
emission is assumed to take place at the redshift of coalescence. We show that sources at 
redshifts z > 0.5 contribute to a truly continuous stochastic background, while sources 
at redshifts 0.2 < z < 0.5 are responsible for a popcorn noise, with duty cycle of 1 and 
0.1 respectively (Fig. p|. The closure density reaches a maximum of ^Ig^ ~ 7.3 x 10~^° 
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Figure 2. closure density of the continuous background produced by DNS coalescences 
at 2; > 0.5 (continuous line) and of the popcorn contribution corresponding to sources 
between z = 0.2 — 0.5 (dashed line). The signal from the whole population down to 
z = is also plotted for comparison (dotted line). 




Figure 3. duty cycle from sources up to redshift z. The horizontal lines at = 1 
and D = 0.1 show the limit of the continuous and the popcorn backgrounds. 
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around 480 Hz for the continuous contribution and of ~ 9.7 x 10~^° around 520 
Hz for the popcorn background (Fig. |2]). The total background, including the nearest 
sources down to 2 ~ is slightly higher, with a maximum of Qg^ ~ 1.2 x 10~^ at 560 Hz. 
^ used a similar procedure to calculate the background in the frequency range of LISA. 
At lower frequencies our results are comparable (with Qg^ ~ 4 — 5 x 10~^° at around 
100 Hz for the total background) besides different assumptions about the SFR, the mass 
range for NS progenitors and the distribution of the coalescence time. But in [9j the 
maximum occurs at lower frequencies (~ 100 Hz), since being interested in the signal in 
the range between 10 /iHz and 1 Hz, they have set the value of the maximum frequency 
to that expected at a separation three times the last stable orbit (z/max = O.lQz/^so). 
Those authors have stressed that their calculations are expected to be accurate in the 
frequency band of LISA, thus a direct comparison with our predictions for the frequency 
band of ground based detectors is probably not very meaningful. 



3. Detection 



The optimal strategy to search for a gaussian (or continuous) stochastic background, 
which can be confounded with the intrinsic noise background of the instrument, is to 
cross correlate the measurements Si of multiple detectors. When the background is 
assumed to be isotropic, unpolarized and stationary, the cross correlation product is 
given by p6|: 

Si*{f)Q{f)s2{f)df (17) 

-oo 

where 



is a filter that maximizes the signal to noise ratio (S/R). In the above equation, 
Pi(/) and P2{f) are the power spectral noise densities of the two detectors and T is the 
non-normalized overlap reduction function, characterizing the loss of sensitivity due to 
the separation and the relative orientation of the detectors (see Fig. |4]). The optimized 
S/N ratio for an integration time T is given by 

■n' 16vr4 Jo ^pPiimifY 

The S/N for the main terrestrial interferometer pairs, at design sensitivity and in 
their advanced configuration, after one year of integration, are given in Table [ij for a 
detection rate a = 90% and a false alarm rate /? = 10%. Expressions for the power 
spectral densities of actual detectors can be found in [31] (see Fig. [s]). 

The continuous signal is below the sensitivity that can be obtained by cross- 
correlating actual pairs of detectors For example, considering co-located and co- 
aligned interferometers, such as Virgo or LIGO, we find a maximum signal-to- noise 



74 .oo r^fflO^ (f) 

/o ^/ fap np;;) M^"^(2/^) - erfc-(2«))-. (19) 
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Figure 4. overlap reduction function for the most promising detector pairs. LHO and 
LLO stand for LIGO Hanford Observatory and LIGO Livingston Observatory. 




Figure 5. designed sensitivities of the main first generation interferometers 
(continuous), compared to the planned sensitivities of the advanced interferometer 
LIGO Ad and the third generation interferometer EGO. 
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initial 


2.8 X 10-3 


5.5 X 10-6 


6.1 X 10-*^ 


4.9 X 10-6 




advanced 


0.52 


0.029 








3^d generation 
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Table 1. Expected signal-to-noise ratio, corresponding to the continuous background 
{D > 1) and for the actual and future terrestrial interferometer pairs for an integration 
time T = 1 year, a detection rate a — 90% and a false alarm rate f3 = 10%. LHO and 
LLO stand for LIGO Hanford Observatory and LIGO Livingston Observatory. 



ratio of S/N ~ 0.003 {S/N ~ 0.5) for the initial (advanced) configuration. However, 
the sensitivity of the future third generation of detectors such as EGO, presently in 
discussion, could be high enough to gain one order of magnitude in the expected signal 
to noise ratio {S/N ~ 8) . On the other hand, the popcorn noise contribution could 
be detected by new data analysis techniques currently under investigation, such as the 
maximum likelihood statistic [29], or methods based on the Probability Event Horizon 
concept [30j, which describes the evolution, as a function of the observation time, of the 
cumulated signal throughout the Universe. 



4. Simulations of the DNS population 

In this section, we introduce Monte Carlo simulations of the extra-galactic population of 
DNSs and calculate the resulting stochastic background. This method can be extended 
to any kind of sources, in particular to GW events that are delayed with respect to the 
formation of the progenitors. The simulations follow the evolution of the system from 
the birth of the progenitors to the merging of the two neutron stars, after the redshift of 
formation and the coalescence time have been selected. The difference with the previous 
simulations of [T^l US] , in addition to the update of the initial mass function, the star 
formation rate and the cosmological model, is that the normalization (the ratio between 
the real number of DNSs and the number of simulated DNSs or runs) is done in a more 
realistic way by considering comoving volume elements instead of redshift intervals when 
following the evolution of the progenitors. 

To simulate a population of coalescences observed today in an element of comoving 
volume, we proceed as follow (see Fig. |6]): 

(i) The time of formation of the progenitors is selected from the probability distribution 

P (T)= ^/^^-^^ (20) 

defined by normalizing in the interval — 12.5 Myr the formation rate of the 
progenitors (eq. [?]) , 

(ii) The cosmic time Tf, at which the progenitors have evolved into two neutron stars 
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Figure 6. flowchart of the Monte-Caiio simulations described in section 3. Each 
system is generated with a cosmic time of formation Tf, (from which the cosmic time 
of formation of the DNS is calculated), and a coalescence time Tc, which defines the 
coalecence cosmic time Tc- Only DNSs which coalesce at redshifts Zc < z^, contribute 
to the integrated signal; their fluences are calculated and combined with adequate 
normalization factors (see text) to compute the total flux and the density parameter 
^gw{i^o)- The critical redshift to have a continuous stochastic background is = 0.4 
(z* = 0.2 for the popcorn noise). 



and start to coalesce is given by 

n = Tf- n (21) 

where Tf, (~ 10® yr) has been defined in the previous section as the mean hfetime 
of the progenitors. 

(iii) The coalescence timescale which depends on both the orbital parameters and the 
masses of the two neutron stars, is selected from the probability distribution eq. [8} 
between 0.2 Myr and 20 Gyr. 

(iv) The cosmic time Tc at which the coalescence occurs is given by 

Tc = Tb- (22) 
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and the corresponding coalescence redshift Zc is derived by solving the equation: 

^'^1, Ho{l + z)E{z) ^^^^ 

(v) Each DNS, thus generated, is then sorted into bins of cosmic time [T^;T^ + ATc], 
for which we calculate the total flux as the sum of all the individual fiuences, 
normalized by the ratio between the total formation rate of the progenitors in the 
range — 12.5 Gyr {N°) and the number of simulated DNSs (Nsim)'- 

N° 

FAn.^o) = ^Y.fiTl.yo) (24) 

' sim j 

with Tl in the range [T^ ; + ATc] and where 

N; = / R%Tf)dTf (25) 

(vi) The model of the star formation rate being isotropic, the element of comoving 
volume at cosmic time T is considered as representative of the entire population 
in the shell between [T; T + dT]. Therefore, the total flux from sources located 
between [T^;T^ + ATc] writes: 

^(^o) = E^.(^c^^o)^(Tc^) (26) 

or equivalent ly 

A^" ■ dV ■ 

F{uo) = ^T.fiT'.^^o)-^iTc) (27) 

-'"'sim j "-^ 

The closure density (eq. [2]) corresponding to the continuous background {z^ > 0.5) 
is plotted in figure (Fig. [?]) and compared with the results obtained in section 1. For 
a number of runs Ngi^ = 10^, the agreement is better than 99.5%, which is accurate 
enough to vahdate the Monte Carlo procedure. 



5. Conclusions 



We presented Monte Carlo numerical simulations of the extra-galactic population 
of double neutron stars and investigated its contribution to the gravitational wave 
stochastic background. The stochastic background formed by the final stage of the 
coalescence in the frequency band 10 — 1500 Hz is continuous for sources beyond 
z ~ 0.5 and rather a popcorn noise between 0.2 < z < 0.5. The closure density of 
the continuous contribution reaches a maximum of Qgw ~ 7.3 x 10^^*^ at around 480 
Hz, which is below the sensitivity of actual and advanced interferometers but may be 
detectable with the third generation. The popcorn contribution seems more promising 
with Qg^ ~ 9.7 X 10~^° at 520 Hz. The advantage of the Monte Carlo simulations, 
compared to the direct calculation are numerous. On the one hand they permit to 
study the statistical properties of the background in the time domain, in particular the 
non gaussian popcorn contribution, for which adapted detection strategies are currently 
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Figure 7. closure density of the continuous background produced by DNS coalescences 
at z > 0.5 derived from the direct calculations described in section 1 (continuous line) 
and from the Monte Carlo simulations described in section 2 (square) , for a number of 
runs Nsim — 10®. The agreement between the two is better than 99.5%. 



under development. The simulation of GW time series [321 [33] that can be injected in 
the output of a pair of detectors |3ll |35] , is essential to test and validate data analysis 
pipelines. On the other hand incorporating new parameters (such as the eccentricity, 
the orbital separation and the masses of the two stars) can be done in a very simple way. 
In this work all the initial informations are included in the coalescence time. However 
to investigate the stochastic background formed by the low frequency inspiral phase, 
in the frequency domain of the spatial detector LISA, when the system can be highly 
eccentric [36l l37] emitting GW at higher harmonics to the keplerian frequency [2T1 [38| |8] , 
one needs to follow the combined evolution of the frequency, the eccentricity and the 
redshift of emission. This work is currently in progress and will be reported in a future 
paper. 
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